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Ten-place tables of spherical Bessel functions of complex argument with real and 
Imaginary values from 0.00] to 1000 In varying Increments for 100 consecutive orders 
fare presented. The method used to cenwrate these tables Is based on a continued 
| fraction representation of the ratio of Bessel functions of consecutive order. This 
method retains a high degree of accuracy In a finite digit machine, Incorporates 
a terminating condition, Is fast and sImle, and Is versatile enough to find wide 
application as a method for continued fractions In general. 
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INTRODUCTION 


Spherical Bessel functions of comp!ex argument find application 

In sclentific computations ranging from electronic circult analysis 

to the understanding of nuclear and Mie scattering. Yet tables of 

Bessel functions are presently |Imited to spherical Bessel functions 

of real argument and to nonspherical Besse! functions of complex eo 
argument of 0 and | order. Tables of the former may be found In 

Abramowltz Cl], and tables of the latter in towan [2]. The 

generation of spherical Bessel functions of complex argument Is not 4 
a trivial problem, even with modern high-speed computers, because Fi 
of difficulties assoclated with each of the methods used for their ' 
generetion. i 


To surmount these difficulties, It became necessury tc develop a 

method which would be simple and accurate for machine computation, { 
and which would Incorporate a terminating condition, The result was q 
a method based on a continued fraction representation of the ratlo of 
Bessel functions of consecutive order. 


The purpose of this report is primarily to present the new method 
of calculating continued fractions, and secondarily to present 
tables of spherical Bessel functions of complex argument against 
which this method and others may be checked. The development 

Is placed In perspective by reviewing the advantages and dis- 
advantages of the other methods and algorithms which might be used, 
Suggestions are made relative to the extension of the method to 
other functions of complex argument. 


RECURSION RELATIONS 


Recursion relations are the simplest and mst common method of 
generating Bessel functions. The fundamental recursion relation Is 


Sq OO) + Jy) = 2v/x J Cx) ap 


where J Is the Bessel function of order v and complex argument x 

and v = m1/2 where n is an Integer. This relation may be used with 
two consecutive orders to generate all the necessary orders by starting 
at low order and recurring upward or starting at an approximate high 
order and recurring downward. 
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Unfortunately both of these methods have serlous [imitatlons. Lowan ; 
[2] has noted that the upward recursions are unstable when 2v/x Is 
greater than unity, because the values of consecutive orders of 
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Bessel functions have common tnitlal digits. Gautch! [3] has shown \ 
that the upward recursion has no digit accuracy after only seven 
recirsions with 32-bit capacity for x = 1.0. Therefore, even using 
double precision arithmetic, the procedure will not work for more 
than a very few consacutive orders. 


Miller's downward recursion scheme [4] has solved the problem of 
accuracy at the expense of machine time and storage. By starting 
with an approximate value for high orders, one can recur downward i 
with only roundoff error. Since downward recursions start at a very 4 
hign order for moderate arguments, this method is impractical 
when only a few orders are needed. This technique also suffers 


from the same problem that the forward recursion has - once - 
accuracy is lost in a term, all the following terms have lost the a 
same amount of accuracy. In this case, roundoff error may accumulate 4 
so that the more significant smaller orders are very much more : 


Inaccurate than ‘the approximate initial high order. 
INFINITE SERIES 


The infinite series for Bessel functions {s se] dom used because the 
serles |s comutationally divergent for moderate x and small orders. 
When the series does converge, it converges slowly, end it is difficult 
to program. Infinite serles represenvations of Bessel functions are 
therefore mathematically Interesting but not very practical for 
computation. 


% 
£ 
# 
“4 


CONTINUED FRACTIONS 


Continued fractions have the advantage of being simple, gasy to program, 
and accurate for machine computation. Unfortunately, they have been 
neglected If favor of recursion relations because 01 difficulttes In 
determining the accuracy of a given number of terms in a continued 
fraction. The following development of continued fractions is !ntended 
to review their basic properties, point out the weaknesses in available a 
algorithms, and to develop a method of generating Bessel functions which Ze 
is free of these weaknesses. 


Basic Properties 


A continued fraction may be written as 


f= a + by 


c, * by 
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= ene {2) 
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which can be rewritten as 


2 ay ee (3) 


The nth approximate convergent to the continued fraction representation 
of f is written as 


pte te (4) 


y Continued fractions may be generated from a three-term recursion 
f relation in a simple way: 
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This form is one of a simple continued fraction which may be defined as 
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where a, # 0 and the a's may be negative. 


Equation (7) may be written in the more convenient notation of 


f(x) = [a,, apr az, See] 


(8) 


The nth convergent is written In | ike manner: 


f(x) = [a,, Ons Bzs vee y aad las 


It has been shown by Wall [51] that all continued fractions may be 
expressed as simple continued fractions of the above form. In 
addition, Watson [6] has shown that the form of the continued fraction 
for Bessel functions is unique, requiring that al! the oth convergents 
of the continued fraction in whatever form have the same numerical 
value. For example, all of the foliowing are equivalent in the nth 
convergent: 
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Algorithms 


it is now possible to generate the ratio of Bessel functions by starting " 
with the last value on the right of Eq. (9) and working backwards 4 
to a,. However, the minimum necessary value of n must be specified uf 
befote the computation starts, since adding one term means the re- . 


computation of the entire fraction. Gautch! [3] has made an estimate 4 
of the value of n for determining Bessel functions to a given accuracy. s 
This method does not apply In general to continued fractions and is i 


only an approximate estimate which is difficult to apply. Moreover, i 
the tail-to-head method suffers somewhat from roundoff error when 
the difference between two values fs near zero. 4 
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Another method of computing Wyay/dy follows logically from the 
relations 
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where p_ is the nin partial numerator and qn the nih partial 
denominator of *he nth convergent. 


‘A Walt [5] has shown that these relations are generally true for any 

{ continued fraction. Starting with the initial values of p, = a,, 

, q) = by Po = 1, and G.. 0, one can generate successive approximants 

g to the exact infinite continued fraction until the desired accuracy Is 
obtainea, This method ts illustrated in Teichroew [7] and Abramowitz 
[1iJ. The former has noted that the values of p. and q. may grow so 
large as to cause overflow in the computer. Moreover, for negative 
aan Initial number of digits may cancel, causing loss of accuracy 
in the mth convergent in exactly the same way as the original 
recursion relation (Eq. (1)). This forward ratio method leaves 
much to be desired, 
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By noticing that the relations in Eq. (13) are themselves recursion 
relations, one obtains a new method of generating Bessel functions: 
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So the value of p,/4, ts given by 
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Eq. (17) gives the value of the nt approximant to the exact continued 
fractions representation for a particular continued fraction. When 
the ratio 
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where € iS the fraximum allowable errcr, the product can be terml nated 
Since al, subsequent ratios will be unity. Eq. (17) gives a simple 
method of computing the ratios of Bessel functions cf successive order 
to the number of significant digits in the computer. The specific 
formula for the Besse! function is 


J - | | 
5 = 2Vx hb RET prea PP sale eee 419) 
Vv -2(vtl)x + 2(vt2)~% 9 + -2(vt3)x + 


Every other a, is negative in Eq. (17), but the application is 
unchanged. Each numerator and denominator is generated from the 
preceding one by taking the reciprocal and adding the next a, 


C2tve2m)x” | ~2(ve2me 1x7) baie avx7 7 rose Covx" !y 
(20) 
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Since both the numerators and denominators alternate in sign, it is 
simpler to consider all a_ positive, and to generate each successive 
term by taking the reciprocal and subtracting from the next a. The 
ath approximate convergent is the product of all the preceding 


ule and the fraction terminates when the nt) ratio In Eq. (18) 
s unity. 


\t may be noted that the order of the terms in the simple continued 
fraction has been reversed between Eqs. (8) and (17). This means that 
it Is not necessary to decide how many terms are desired before the 
computation, as tn Eq. (9). Because the successive numerators and 
denominators of Eq. (17) are tending toward the same finite value, 
there is no large growth in the numerator and denominator as in the 
use of Eq. (13). 


An advantage that is not Immediately apparent Is that a particular 

n convergent being negr zero or very large is no disadvantage. 
Although the computes n +,epproximant will deviate slightly from the 
theoretical vaiue, the n convergent, wil | converge back to the 
correct value. {In other words, the n approximant Is slightly self- 
correct ing. 


The ccrtinued fraction converges very quickly and does not require 
significant storege. This first ratio usually requires the largest 
number of terms, and for x = 1.0 and v = 1.5 with 64-bit capacity, 
only nine terms are needed. A sample calculation follows: 
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argument tabulated in Abramowitz [i]. The agreement between the 
doubie precision roviine and the handbook was exact for 10 significant 
digits for all arguments. By comparing the single precIsion 

version to the double precision version, it Is astimated that the 

last tour of the 16 digits would be suscten*'..|e@ te facnine roundoff. 
Therefore, the table is reproduced to 10 "o ‘slow for the 
additional possibility of roundoff err “ vl. « arguments. 


Ae 
i x= 1.0 v= 9.5 

| 

= 

Z 

7 el =. (19) (20.94736842) (22.95226131) (24.95643131) (26. 35993017) 
an Jy (21) (22.95238095) (24.95643154) (26,95993017) 
' 

e This value ts identical to the handbook vatue [17]. 

wt 

\ The development has made no restrictions on the form of the argument 
ie or order, so they can both be comlex and positive or negative. 

: Near zeros of the Bessel functions It should be noted that the 

ce accuracy may drop to ten places rather than ten significant digits. 

. The actual spherical Bessel functions rather than their ratios as 

i In Eq. (19) are found from 
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ACCURACY 

7 A conparison between the accuracy of a double precision 64-bit and 

‘ single precision 32-bit program was made for the real values of 
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a OTHER FUNCTICHS 
Many known functions have been expres=qid .2 4. .*rt, of continued 


fractions, and they may be evaluated by Eq, (17). Abramowitz [1] has 
listed one form of a-continued traction re; resentation for mst of the 
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functions In his handbook. The form of any continued fraction may 
be changed to inat of Eq. (7) by the following transformations: 


see (23) 


The only requirement is that a_ be nonzero; if a_ ts zero, this 


requirement may be met by takIRg the reciprocal 8f the continued 
fraction. 


The continued fraction technique will converge with fewer terms 

than any other current method and Is specifically adapted to machine 
computation. Although mathematically equivalent to recursion 
relations, the continued fraction technique gains accuracy from the 
fact that the product of two numbers In the computer does not lose 
accuracy, as does the difference of two numbers. 
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APPENDIX A 
USE OF THE TABLES 


The fo!lowing tables contain [00 orders of spherical Bessel functions 
of the first and second «ind for the same argument on the same page. 
The real and Imaginary parts of the spherical Bessel functions of 

the first kind f(x) = Van] 416%) are presented on the left side 


with the appresrfate order_n, and the spherica! Bessel functions of 
the second king y, (x) = Pam XY ga (x)= (= 19" Sena (x) are in 


-(nt4g) 
like manner on the right side of the page. 

Beyond the 20th order, the missing orders can be found by a down- 
ward recursion for /_(z) and by an upward recursion for y (x). 

The spherical Besse!" functions of the complex conjugate of the 
argument presented here can be found by raking the comlex conjugate 
of the results in the Table. Since these Tables are intended as 


check values, Interpolation should not be attempted between values 
of argument. 


In those fs cases where the value 98 the spherical Bessel functions 
exceed 107-00 g scale factor of 10490 was multiplied by the fi (x) 
and 107200 was multiplied by the y (x). These scale factors are 
indicated by a sudden break in the order of magnitude with increasing 
n since the initial values never exceeded the scale factor. 


Many of the more useful relations among the spherical Bessel functions 
may be found in Abramowitz [I]. If these relations are used to, 

check the table, It should be noted that the Initial values of J 

and yare given by Eq. (22) and that the last or highest order being 
correct depends on all the lower orders being correct. 
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